##
## Government and Opposition in Legislative Speechmaking: 
## Using Text-As-Data to Estimate Brazilian Political Parties Policy Positions
##
## OPINION IRT MODEL (R-32 bits)
##

#working directory
setwd("C:\\Users\\mauricio\\Desktop\\opinion_irt\\")

#load packages
require(tm)
require(SnowballC)
require(e1071) #naive bayes
require(sqldf)
require(rjags)
require(ggplot2)

#function to remove leading and trailing whitespace
trim <- function( x ) {
  gsub("(^[[:space:]]+|[[:space:]]+$)", "", x)
}

## Function to concatenate lines in a list
concat <- function(v){
	res = ""
	for(i in 1:length(v))
		res <- paste(res, v[i], sep = " ")
		res
}

#############################################################

##
## 01. Sentiment analysis (Naive bayes)
##

#load data and merge
disc <- read.csv("disc_95-14.csv")

training <- readRDS("training_set.rds")
training <- training[,c("link", "sentiment")]

disc <- merge(disc, training, all.x = T, by = "link")

##
## cleaning resumo
##

##resumo (factor to character)
disc$resumo <- as.character(disc$resumo)

##Lowercase
disc$resumo <- tolower(disc$resumo)

##Remove punctuation
disc$resumo <- removePunctuation(disc$resumo)

##Remove numbers
disc$resumo <- removeNumbers(disc$resumo)

##Remove stopwords
disc$resumo <- removeWords(disc$resumo, stopwords(kind = "portuguese"))

##Stemming
pron <- list(NULL)
for(i in 1:length(disc$resumo)){
	pron[[i]] <- wordStem(strsplit(disc$resumo[i], " ")[[1]], language = "portuguese")
	pron[[i]] <- concat(pron[[i]])
}

teste <- unlist(pron)
disc$resumo <- teste

##Eliminating extra space
disc$resumo <- stripWhitespace(disc$resumo)

###

##Create Corpus
corpus <- VCorpus(VectorSource(disc$resumo))

##Document term matrix
ndocs <- length(corpus)
minDocFreq <- ndocs * 0.01 #ignore sparse terms (apperaing in less than 1% of the docs)
maxDocFreq <- ndocs * 0.99 # ignore common terms (appearing in more than 99% of the docs)
dtm <- DocumentTermMatrix(corpus, control = list(bounds = list(global = c(minDocFreq, maxDocFreq))))
dtm <- as.matrix(dtm)

dtm <- replace(dtm, dtm > 1, 1)

###

## Train the model
training_set <- which(is.na(disc$sentiment) == F)
classifier <- naiveBayes(dtm[training_set,], as.factor(disc$sentiment[training_set]))

## Test the model
test_set <- which(is.na(disc$sentiment) == T)
predicted <- predict(classifier, dtm[test_set,])

teste <- disc[test_set,]
teste$sentiment <- predicted

##Accuracy e evaluation

#DO NOT RUN AGAIN!
#IT CREATED THE TABLE FOR ACCURACY EVALUATION
#AFTER IT, WE HAND-CODED THE SPEECHES FOR CROSS-VALIDATION
#acc <- sample(c(1:nrow(teste)), 1000)
#accuracy <- teste[acc, c(1:5,8:10)]
#
#temp <- read.csv("disc_95-14.csv")[,c("link", "resumo")]
#accuracy <- merge(accuracy, temp, all.x = T, by = "link")
#
#write.csv(accuracy, "accuracy.csv", row.names = F)

##
## TABLE 3
##

accuracy <- readRDS("accuracy.rds")
table(accuracy$handcode, accuracy$sentiment)
(384+378)/10
#Accuracy of 76.2%

#####################################################################

##
## 2. Statistical model
##

training <- read.csv("training_set.csv", sep = ";")
training <- training[,c("link","senador","partido","uf","data",
				"catalogo","ano","sentiment")]
teste <- teste[,c("link","senador","partido","uf","data",
				"catalogo","ano","sentiment")]
dados <- rbind(training, teste)

#Recoding partido
dados$partido <- ifelse(dados$partido %in% c("PL", "PR"), "PL>PR",
			as.character(dados$partido))
dados$partido <- ifelse(dados$partido %in% c("PPB", "PPR", "PP"), "PDS>PP",
			as.character(dados$partido))
dados$partido <- ifelse(dados$partido %in% c("PMR", "PRB"), "PMR>PRB",
			as.character(dados$partido))
dados$partido <- ifelse(dados$partido == "DEM", "PFL>DEM",
			as.character(dados$partido))

##
## Dilma 1
##

#Subset data
tab <- subset(dados, dados$ano >= 2011)

#Subset parties (at least 3 senators and 500 speeches)
table(tab$partido)

for(i in unique(tab$partido)){
	print(i)
	print(length(unique(subset(tab$senador, tab$partido == i))))
}

#tab <- subset(tab, tab$partido %in% c("PSDB", "PMDB", "PDS>PP", "PT", "PSB", "PTB", "PDT"))
tab <- subset(tab, tab$partido %in% c("PSDB", "PMDB", "PDS>PP", "PT", "PSB", "PTB",
						"PDT", "PFL>DEM", "PL>PR"))

#categories
temp <- unlist(strsplit(as.character(tab$catalogo), "[.]"))
temp <- trim(temp)
temp <- subset(temp, temp != "")
temp <- subset(temp, temp != "D")
temp <- subset(temp, temp != "P")
temp <- gsub(" ", "", temp)
temp <- paste("\\b", temp, "\\b", sep = "")
categ <- temp
sort(unique(categ))

#merge index in categories
tab$catalogo <- gsub(" ", "", tab$catalogo)
tab$catalogo <- gsub("[.]", " ", tab$catalogo)
tab$catalogo <- trim(tab$catalogo)

#merge
tab <- tab[,c("senador", "partido", "catalogo","sentiment")]

temp <- list(NULL)
categoria <- list(NULL)
k <- 0
for(i in unique(categ)){
	k <- k + 1
	temp[[k]] <- tab[grep(i, tab$catalogo),]
	categoria[[k]] <- rep(i, nrow(temp[[k]]))
}

final <- do.call("rbind", temp)
final$categoria <- unlist(categoria)
final$categoria <- gsub("[\\]b", "", final$categoria)

##
## group data by party
##

final <- sqldf("SELECT partido, categoria, sum(sentiment) AS npos, Count(*) AS tot
		FROM final
		GROUP BY partido, categoria")

##
## create matrix with "n" and "p" (binomial stochastic component)
##

s <- reshape(final[,c(1,2,3)],
	v.names = "npos",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn <- reshape(final[,c(1,2,4)],
	v.names = "tot",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn[is.na(nn)] <- 10

#polarities
s <- rbind(s[which(s$partido == "PT"),], s[-which(s$partido == "PT"),]) #pt
nn <- rbind(nn[which(nn$partido == "PT"),], nn[-which(nn$partido == "PT"),]) #pt
s <- rbind(s[which(s$partido == "PSDB"),], s[-which(s$partido == "PSDB"),]) #psdb
nn <- rbind(nn[which(nn$partido == "PSDB"),], nn[-which(nn$partido == "PSDB"),]) #psdb

#######################################
#######################################

##
## JAGS MODEL
##

N <- nrow(s)
M <- ncol(s[,-1])
y <- s[,-1]
n <- nn[,-1]

forJags <- list("y" = y, 
	"n" = n,
	"N" = N, 
	"M" = M)

jagsModel <- jags.model(file = "binomial_model.bug",
	data = forJags)

binomModel_dilma <- coda.samples(jagsModel,
	variable.names = c("alpha", "beta", "theta"),
	n.iter = 50000,
	burnin = 5000,
	thin = 50)

##
## Summary
##

summary(binomModel_dilma)

##################################################################################
##################################################################################
##################################################################################
##################################################################################

##
## Lula 2
##

#Subset data
tab <- subset(dados, dados$ano >= 2007 & dados$ano <= 2010)
table(tab$partido)

#Subset parties (at least 3 senators)
for(i in unique(tab$partido)){
	print(i)
	print(length(unique(subset(tab$senador, tab$partido == i))))
}

#tab <- subset(tab, tab$partido %in% c("PSDB", "PT", "PMDB", "PFL>DEM", "PDT", "PTB", "PL>PR"))
tab <- subset(tab, tab$partido %in% c("PMR>PRB", "PSDB", "PT", "PMDB", "PSB",
						"PFL>DEM", "PDT", "PTB", "PL>PR"))

#categories
temp <- unlist(strsplit(as.character(tab$catalogo), "[.]"))
temp <- trim(temp)
temp <- subset(temp, temp != "")
temp <- subset(temp, temp != "D")
temp <- subset(temp, temp != "P")
temp <- gsub(" ", "", temp)
temp <- paste("\\b", temp, "\\b", sep = "")
categ <- temp
sort(unique(categ))

#merge index in categories
tab$catalogo <- gsub(" ", "", tab$catalogo)
tab$catalogo <- gsub("[.]", " ", tab$catalogo)
tab$catalogo <- trim(tab$catalogo)

#merge
tab <- tab[,c("senador", "partido", "catalogo","sentiment")]

temp <- list(NULL)
categoria <- list(NULL)
k <- 0
for(i in unique(categ)){
	k <- k + 1
	temp[[k]] <- tab[grep(i, tab$catalogo),]
	categoria[[k]] <- rep(i, nrow(temp[[k]]))
}

final <- do.call("rbind", temp)
final$categoria <- unlist(categoria)
final$categoria <- gsub("[\\]b", "", final$categoria)

##
## group data by party
##

final <- sqldf("SELECT partido, categoria, sum(sentiment) AS npos, Count(*) AS tot
		FROM final
		GROUP BY partido, categoria")

##
## create matrix with "n" and "p" (binomial stochastic component)
##

s <- reshape(final[,c(1,2,3)],
	v.names = "npos",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn <- reshape(final[,c(1,2,4)],
	v.names = "tot",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn[is.na(nn)] <- 10

#polarities
s <- rbind(s[which(s$partido == "PT"),], s[-which(s$partido == "PT"),]) #pt
nn <- rbind(nn[which(nn$partido == "PT"),], nn[-which(nn$partido == "PT"),]) #pt
s <- rbind(s[which(s$partido == "PSDB"),], s[-which(s$partido == "PSDB"),]) #psdb
nn <- rbind(nn[which(nn$partido == "PSDB"),], nn[-which(nn$partido == "PSDB"),]) #psdb

#######################################
#######################################

##
## JAGS MODEL
##

N <- nrow(s)
M <- ncol(s[,-1])
y <- s[,-1]
n <- nn[,-1]

forJags <- list("y" = y, 
	"n" = n,
	"N" = N, 
	"M" = M)

jagsModel <- jags.model(file = "binomial_model.bug",
	data = forJags)

binomModel_lula2 <- coda.samples(jagsModel,
	variable.names = c("alpha", "beta", "theta"),
	n.iter = 50000,
	burnin = 5000,
	thin = 50)

##
## Summary
##

summary(binomModel_lula2)

##################################################################################
##################################################################################
##################################################################################
##################################################################################

##
## Lula 1
##

#Subset data
tab <- subset(dados, dados$ano >= 2003 & dados$ano <= 2006)
table(tab$partido)

#Subset parties (at least 3 senators)
for(i in unique(tab$partido)){
	print(i)
	print(length(unique(subset(tab$senador, tab$partido == i))))
}

#tab <- subset(tab, tab$partido %in% c("PSDB", "PMDB",	"PFL>DEM", "PT", "PDT"))
tab <- subset(tab, tab$partido %in% c("PSDB", "PSB", "PTB", "PMDB",
						"PFL>DEM", "PPS", "PT", "PDT", "PL>PR"))

#categories
temp <- unlist(strsplit(as.character(tab$catalogo), "[.]"))
temp <- trim(temp)
temp <- subset(temp, temp != "")
temp <- subset(temp, temp != "D")
temp <- subset(temp, temp != "P")
temp <- gsub(" ", "", temp)
temp <- paste("\\b", temp, "\\b", sep = "")
categ <- temp
sort(unique(categ))

#merge index in categories
tab$catalogo <- gsub(" ", "", tab$catalogo)
tab$catalogo <- gsub("[.]", " ", tab$catalogo)
tab$catalogo <- trim(tab$catalogo)

#merge
tab <- tab[,c("senador", "partido", "catalogo","sentiment")]

temp <- list(NULL)
categoria <- list(NULL)
k <- 0
for(i in unique(categ)){
	k <- k + 1
	temp[[k]] <- tab[grep(i, tab$catalogo),]
	categoria[[k]] <- rep(i, nrow(temp[[k]]))
}

final <- do.call("rbind", temp)
final$categoria <- unlist(categoria)
final$categoria <- gsub("[\\]b", "", final$categoria)

##
## group data by party
##

final <- sqldf("SELECT partido, categoria, sum(sentiment) AS npos, Count(*) AS tot
		FROM final
		GROUP BY partido, categoria")

##
## create matrix with "n" and "p" (binomial stochastic component)
##

s <- reshape(final[,c(1,2,3)],
	v.names = "npos",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn <- reshape(final[,c(1,2,4)],
	v.names = "tot",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn[is.na(nn)] <- 10

#polarities
s <- rbind(s[which(s$partido == "PT"),], s[-which(s$partido == "PT"),]) #pt
nn <- rbind(nn[which(nn$partido == "PT"),], nn[-which(nn$partido == "PT"),]) #pt
s <- rbind(s[which(s$partido == "PSDB"),], s[-which(s$partido == "PSDB"),]) #psdb
nn <- rbind(nn[which(nn$partido == "PSDB"),], nn[-which(nn$partido == "PSDB"),]) #psdb

#######################################
#######################################

##
## JAGS MODEL
##

N <- nrow(s)
M <- ncol(s[,-1])
y <- s[,-1]
n <- nn[,-1]

forJags <- list("y" = y, 
	"n" = n,
	"N" = N, 
	"M" = M)

jagsModel <- jags.model(file = "binomial_model.bug",
	data = forJags)

binomModel_lula1 <- coda.samples(jagsModel,
	variable.names = c("alpha", "beta", "theta"),
	n.iter = 50000,
	burnin = 5000,
	thin = 50)

##
## Summary
##

summary(binomModel_lula1)

##################################################################################
##################################################################################
##################################################################################
##################################################################################

##
## FHC 2
##

#Subset data
tab <- subset(dados, dados$ano >= 1999 & dados$ano <= 2002)
table(tab$partido)

#Subset parties (at least 3 senators)
for(i in unique(tab$partido)){
	print(i)
	print(length(unique(subset(tab$senador, tab$partido == i))))
}

#tab <- subset(tab, tab$partido %in% c("PMDB", "PT", "PSDB", "PFL>DEM", "PSB"))
tab <- subset(tab, tab$partido %in% c("PMDB", "PT", "PDT", "PSDB", "PDS>PP",
						"PFL>DEM", "PPS", "PSB", "PTB"))

#categories
temp <- unlist(strsplit(as.character(tab$catalogo), "[.]"))
temp <- trim(temp)
temp <- subset(temp, temp != "")
temp <- subset(temp, temp != "D")
temp <- subset(temp, temp != "P")
temp <- gsub(" ", "", temp)
temp <- paste("\\b", temp, "\\b", sep = "")
categ <- temp
sort(unique(categ))

#merge index in categories
tab$catalogo <- gsub(" ", "", tab$catalogo)
tab$catalogo <- gsub("[.]", " ", tab$catalogo)
tab$catalogo <- trim(tab$catalogo)

#merge
tab <- tab[,c("senador", "partido", "catalogo","sentiment")]

temp <- list(NULL)
categoria <- list(NULL)
k <- 0
for(i in unique(categ)){
	k <- k + 1
	temp[[k]] <- tab[grep(i, tab$catalogo),]
	categoria[[k]] <- rep(i, nrow(temp[[k]]))
}

final <- do.call("rbind", temp)
final$categoria <- unlist(categoria)
final$categoria <- gsub("[\\]b", "", final$categoria)

##
## group data by party
##

final <- sqldf("SELECT partido, categoria, sum(sentiment) AS npos, Count(*) AS tot
		FROM final
		GROUP BY partido, categoria")

##
## create matrix with "n" and "p" (binomial stochastic component)
##

s <- reshape(final[,c(1,2,3)],
	v.names = "npos",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn <- reshape(final[,c(1,2,4)],
	v.names = "tot",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn[is.na(nn)] <- 10

#polarities
s <- rbind(s[which(s$partido == "PT"),], s[-which(s$partido == "PT"),]) #pt
nn <- rbind(nn[which(nn$partido == "PT"),], nn[-which(nn$partido == "PT"),]) #pt
s <- rbind(s[which(s$partido == "PSDB"),], s[-which(s$partido == "PSDB"),]) #psdb
nn <- rbind(nn[which(nn$partido == "PSDB"),], nn[-which(nn$partido == "PSDB"),]) #psdb

#######################################
#######################################

##
## JAGS MODEL
##

N <- nrow(s)
M <- ncol(s[,-1])
y <- s[,-1]
n <- nn[,-1]

forJags <- list("y" = y, 
	"n" = n,
	"N" = N, 
	"M" = M)

jagsModel <- jags.model(file = "binomial_model.bug",
	data = forJags)

binomModel_fhc2 <- coda.samples(jagsModel,
	variable.names = c("alpha", "beta", "theta"),
	n.iter = 50000,
	burnin = 5000,
	thin = 50)

##
## Summary
##

summary(binomModel_fhc2)

##################################################################################
##################################################################################
##################################################################################
##################################################################################

##
## FHC 1
##

#Subset data
tab <- subset(dados, dados$ano >= 1995 & dados$ano <= 1998)
table(tab$partido)

#Subset parties (at least 3 senators)
for(i in unique(tab$partido)){
	print(i)
	print(length(unique(subset(tab$senador, tab$partido == i))))
}

#tab <- subset(tab, tab$partido %in% c("PT", "PMDB", "PSDB", "PFL>DEM", "PDS>PP", "PDT"))
tab <- subset(tab, tab$partido %in% c("PT", "PMDB", "PSDB", "PFL>DEM", "PDS>PP",
						"PTB", "PDT"))

#categories
temp <- unlist(strsplit(as.character(tab$catalogo), "[.]"))
temp <- trim(temp)
temp <- subset(temp, temp != "")
temp <- subset(temp, temp != "D")
temp <- subset(temp, temp != "P")
temp <- gsub(" ", "", temp)
temp <- paste("\\b", temp, "\\b", sep = "")
categ <- temp
sort(unique(categ))

#merge index in categories
tab$catalogo <- gsub(" ", "", tab$catalogo)
tab$catalogo <- gsub("[.]", " ", tab$catalogo)
tab$catalogo <- trim(tab$catalogo)

#merge
tab <- tab[,c("senador", "partido", "catalogo","sentiment")]

temp <- list(NULL)
categoria <- list(NULL)
k <- 0
for(i in unique(categ)){
	k <- k + 1
	temp[[k]] <- tab[grep(i, tab$catalogo),]
	categoria[[k]] <- rep(i, nrow(temp[[k]]))
}

final <- do.call("rbind", temp)
final$categoria <- unlist(categoria)
final$categoria <- gsub("[\\]b", "", final$categoria)

##
## group data by party
##

final <- sqldf("SELECT partido, categoria, sum(sentiment) AS npos, Count(*) AS tot
		FROM final
		GROUP BY partido, categoria")

##
## create matrix with "n" and "p" (binomial stochastic component)
##

s <- reshape(final[,c(1,2,3)],
	v.names = "npos",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn <- reshape(final[,c(1,2,4)],
	v.names = "tot",
	idvar = "partido",
	timevar = "categoria",
	direction = "wide")

nn[is.na(nn)] <- 10

#polarities
s <- rbind(s[which(s$partido == "PT"),], s[-which(s$partido == "PT"),]) #pt
nn <- rbind(nn[which(nn$partido == "PT"),], nn[-which(nn$partido == "PT"),]) #pt
s <- rbind(s[which(s$partido == "PSDB"),], s[-which(s$partido == "PSDB"),]) #psdb
nn <- rbind(nn[which(nn$partido == "PSDB"),], nn[-which(nn$partido == "PSDB"),]) #psdb

#######################################
#######################################

##
## JAGS MODEL
##

N <- nrow(s)
M <- ncol(s[,-1])
y <- s[,-1]
n <- nn[,-1]

forJags <- list("y" = y, 
	"n" = n,
	"N" = N, 
	"M" = M)

jagsModel <- jags.model(file = "binomial_model.bug",
	data = forJags)

binomModel_fhc1 <- coda.samples(jagsModel,
	variable.names = c("alpha", "beta", "theta"),
	n.iter = 50000,
	burnin = 5000,
	thin = 50)

##
## Summary
##

summary(binomModel_fhc1)

#################################
#################################

##
## Summaries (mean, lower, upper, sd)
##

##
## FHC 1
##

alpha <- binomModel_fhc1[[1]][,grep("^alpha", dimnames(binomModel_fhc1[[1]])[[2]])]
beta <- binomModel_fhc1[[1]][,grep("^beta", dimnames(binomModel_fhc1[[1]])[[2]])]

theta <- binomModel_fhc1[[1]][,grep("^theta", dimnames(binomModel_fhc1[[1]])[[2]])]
idealpt <- apply(theta, 2, mean)
q.theta <- apply(theta, 2, function(x) quantile(x, probs = c(.025, .975)))
lower <- q.theta[1,]
upper <- q.theta[2,]
sd.theta <- apply(theta, 2, sd)
partido <- c("PSDB", "PT", "PDS>PP", "PDT", "PFL>DEM", "PMDB", "PTB")

tab.fhc1 <- data.frame(partido, idealpt, lower, upper, sd.theta)

tab.fhc1 <- tab.fhc1[order(tab.fhc1$idealpt),]
tab.fhc1$y <- c(1:nrow(tab.fhc1))
tab.fhc1$term <- "Cardoso 1st term"

##
## FHC 2
##

alpha <- binomModel_fhc2[[1]][,grep("^alpha", dimnames(binomModel_fhc2[[1]])[[2]])]
beta <- binomModel_fhc2[[1]][,grep("^beta", dimnames(binomModel_fhc2[[1]])[[2]])]

theta <- binomModel_fhc2[[1]][,grep("^theta", dimnames(binomModel_fhc2[[1]])[[2]])]
idealpt <- apply(theta, 2, mean)
q.theta <- apply(theta, 2, function(x) quantile(x, probs = c(.025, .975)))
lower <- q.theta[1,]
upper <- q.theta[2,]
sd.theta <- apply(theta, 2, sd)
partido <- c("PSDB", "PT", "PDS>PP", "PDT", "PFL>DEM", "PMDB", "PPS", "PSB", "PTB")

tab.fhc2 <- data.frame(partido, idealpt, lower, upper, sd.theta)

tab.fhc2 <- tab.fhc2[order(tab.fhc2$idealpt),]
tab.fhc2$y <- c(1:nrow(tab.fhc2))
tab.fhc2$term <- "Cardoso 2nd term"

##
## Lula 1
##

alpha <- binomModel_lula1[[1]][,grep("^alpha", dimnames(binomModel_lula1[[1]])[[2]])]
beta <- binomModel_lula1[[1]][,grep("^beta", dimnames(binomModel_lula1[[1]])[[2]])]

theta <- binomModel_lula1[[1]][,grep("^theta", dimnames(binomModel_lula1[[1]])[[2]])]
idealpt <- apply(theta, 2, mean)
q.theta <- apply(theta, 2, function(x) quantile(x, probs = c(.025, .975)))
lower <- q.theta[1,]
upper <- q.theta[2,]
sd.theta <- apply(theta, 2, sd)
partido <- c("PSDB", "PT", "PDT", "PFL>DEM", "PL>PR", "PMDB", "PPS", "PSB", "PTB")

tab.lula1 <- data.frame(partido, idealpt, lower, upper, sd.theta)

tab.lula1 <- tab.lula1[order(tab.lula1$idealpt),]
tab.lula1$y <- c(1:nrow(tab.lula1))
tab.lula1$term <- "Lula 1st term"

##
## Lula 2
##

alpha <- binomModel_lula2[[1]][,grep("^alpha", dimnames(binomModel_lula2[[1]])[[2]])]
beta <- binomModel_lula2[[1]][,grep("^beta", dimnames(binomModel_lula2[[1]])[[2]])]

theta <- binomModel_lula2[[1]][,grep("^theta", dimnames(binomModel_lula2[[1]])[[2]])]
idealpt <- apply(theta, 2, mean)
q.theta <- apply(theta, 2, function(x) quantile(x, probs = c(.025, .975)))
lower <- q.theta[1,]
upper <- q.theta[2,]
sd.theta <- apply(theta, 2, sd)
partido <- c("PSDB", "PT", "PDT", "PFL>DEM", "PL>PR", "PMDB", "PMR>PRB", "PSB", "PTB")

tab.lula2 <- data.frame(partido, idealpt, lower, upper, sd.theta)

tab.lula2 <- tab.lula2[order(tab.lula2$idealpt),]
tab.lula2$y <- c(1:nrow(tab.lula2))
tab.lula2$term <- "Lula 2nd term"

##
## Dilma 1
##

theta <- binomModel_dilma[[1]][,grep("^theta", dimnames(binomModel_dilma[[1]])[[2]])]
idealpt <- apply(theta, 2, mean)
q.theta <- apply(theta, 2, function(x) quantile(x, probs = c(.025, .975)))
lower <- q.theta[1,]
upper <- q.theta[2,]
sd.theta <- apply(theta, 2, sd)
partido <- c("PSDB", "PT", "PDS>PP", "PDT", "PFL>DEM", "PL>PR", "PMDB", "PSB", "PTB")

tab.dilma1 <- data.frame(partido, idealpt, lower, upper, sd.theta)

tab.dilma1 <- tab.dilma1[order(tab.dilma1$idealpt),]
tab.dilma1$y <- c(1:nrow(tab.dilma1))
tab.dilma1$term <- "Dilma 1st term"


##############

#RESULTS ("sentiment_irt.csv")
final <- rbind(tab.fhc1, tab.fhc2, tab.lula1, tab.lula2, tab.dilma1)

##################################################

##
## Geweke diagnostic
##

geweke <- as.data.frame(geweke.diag(binomModel_fhc1)[[1]][1]$z)
table(abs(geweke) > 1.96)
geweke <- as.data.frame(geweke.diag(binomModel_fhc2)[[1]][1]$z)
table(abs(geweke) > 1.96)
geweke <- as.data.frame(geweke.diag(binomModel_lula1)[[1]][1]$z)
table(abs(geweke) > 1.96)
geweke <- as.data.frame(geweke.diag(binomModel_lula2)[[1]][1]$z)
table(abs(geweke) > 1.96)
geweke <- as.data.frame(geweke.diag(binomModel_dilma)[[1]][1]$z)
table(abs(geweke) > 1.96)

